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A  tightly  coupled  non-equilibrium  magneto-hydrodynamic  model  for 
Inductively  Coupled  RF  Plasmas 

A.  Munafd,1,a)  S.  A.  Alfuhaid,1, b)  J.-L.  Cambier,2,0)  and  M.  Panesi1,d) 

b Department  of  Aerospace  Engineering,  University  of  Illinois  at  Urbana- Champaign,  Talbot  Lab.,  104  S.  Wright  St., 
Urbana,  61801  IL 

2) Edwards  Air  Force  Base  Research  Lab.,  10  E.  Saturn  Blvd.,  CA  93524 
(Dated:  7  May  2015) 

The  objective  of  the  present  work  is  the  development  a  tightly  coupled  magneto-hydrodynamic  model  for 
Inductively  Coupled  Radio-Frequency  (RF)  Plasmas.  Non  Local  Thermodynamic  Equilibrium  (NLTE)  effects 
are  described  based  on  a  hybrid  State-to-State  (StS)  approach.  A  multi-temperature  formulation  is  used 
to  account  for  thermal  non-equilibrium  between  translation  of  heavy-particles  and  vibration  of  molecules. 
Excited  electronic  states  of  atoms  are  instead  treated  as  separate  pseudo-species,  allowing  for  non-Boltzmann 
distributions  of  their  populations.  Free-electrons  are  assumed  Maxwellian  at  their  own  temperature.  The 
governing  equations  for  the  electro-magnetic  field  and  the  gas  properties  ( e.g .  chemical  composition  and 
temperatures)  are  written  as  a  coupled  system  of  time-dependent  conservation  laws.  Steady-state  solutions 
are  obtained  by  means  of  an  implicit  Finite  Volume  method.  The  results  obtained  in  both  LTE  and  NLTE 
conditions  over  a  broad  spectrum  of  operating  conditions  demonstrate  the  robustness  of  the  proposed  coupled 
numerical  method.  The  analysis  of  chemical  composition  and  temperature  distributions  along  the  torch  radius 
shows  that:  (i)  the  use  of  the  LTE  assumption  may  lead  to  an  inaccurate  prediction  of  the  thermo-chemical 
state  of  the  gas,  and  (ii)  non-equilibrium  phenomena  play  a  significant  role  close  the  walls,  due  to  the  combined 
effects  of  Ohmic  heating  and  macroscopic  gradients. 


I.  INTRODUCTION 

Inductively  Coupled  Plasma  (ICP)  torches  have  wide 
range  of  possible  applications  which  include  deposition 
of  metal  coatings,  synthesis  of  ultra-fine  powders,  gener¬ 
ation  of  high  purity  silicon  and  testing  of  thermal  pro¬ 
tection  materials  for  atmospheric  entry  vehicles.1,2  In  its 
simplest  configuration,  an  ICP  torch  consists  of  a  quartz 
tube  surrounded  by  an  inductor  coil  made  of  a  series  of 
parallel  current-carrying  rings  (see  Fig.  1). 


FIG.  1.  Example  of  ICP  torch  in  operating  conditions  (mini¬ 
torch  facility;  credits  von  Karman  Institute  for  Fluid  Dynam¬ 
ics). 


The  radio- frequency  (RF)  currents  running  through 
the  inductor  induce  toroidal  currents  in  the  gas  which 
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is  heated  thanks  to  Ohmic  dissipation.2,3  If  the  energy 
supplied  is  large  enough,  the  gas  flowing  through  the 
torch  can  undergo  ionization,  leading  to  the  formation 
of  a  plasma. 

The  physico-chemical  modeling  of  the  flow-field  and 
electromagnetic  phenomena  inside  an  ICP  torch  requires, 
in  theory,  the  coupled  solution  of  the  Navier-Stokes  and 
the  Maxwell  equations.  The  numerical  solution  of  this 
coupled  system  of  partial  differential  equations  represents 
a  challenging  task,  due  to  the  disparity  between  the  flow 
and  the  electro-magnetic  field  time-scales.4  Since  in  the 
modeling  of  ICP  facilities  one  is  not  normally  interested 
in  resolving  electro-magnetic  field  oscillations,5  displace¬ 
ment  currents  can  be  safely  neglected  without  introduc¬ 
ing  an  appreciable  error.2,6  This  leads  to  a  more  tractable 
formulation,  as  it  eliminates  the  speed  of  light  from  the 
eigenvalues  of  the  governing  equations.4 

The  first  attempts  to  model  the  temperature  and 
electro-magnetic  field  distributions  inside  ICP  torches 
were  published  in  the  1960-1970’s.  Examples  are  the 
works  of  Freeman  and  Chase,7  Keefer  et.  al.,8  and  the 
series  of  papers  by  Eckert. 9-12  In  most  of  these  references, 
the  torch  was  approximated  as  an  infinite  solenoid  and 
the  plasma  generated  was  considered  in  Local  Thermo¬ 
dynamic  Equilibrium  (LTE)  conditions.  This  reduces  the 
problem  to  the  coupled  solution  of  the  energy  equation 
for  the  gas  (known  as  the  Elenbaas-Heller  equation7,13) 
and  an  induction  equation  for  the  electric  field.  The 
induction  equation  is  formally  identical  to  the  one  de¬ 
scribing  induction  heating  of  metals.14,15  The  authors 
use  the  heat  conduction  potential  (i.e.  s  =  f  ndT),  in 
place  of  the  temperature  T ,  as  thermodynamic  variable. 
This  choice  allows  one  to  hide  the  non-linearity  of  the 
gas  (total)  thermal  conductivity  n  and  can  partially  alle- 
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viate  numerical  instabilities  that  may  arise  when  solv¬ 
ing  the  discretized  set  of  equations  by  means  of  an 
iterative  procedure.  As  recognized  by  Pridmore  and 
Brown,16  the  use  of  the  heat  conduction  potential  be¬ 
comes  less  effective  for  molecular  gases  (e.g.  air  and 
nitrogen  mixtures).  The  major  developments  achieved 
in  the  field  of  Computational  Fluid  Dynamics  (CFD) 
during  the  1970-1980’s,  led  to  the  possibility  of  solv¬ 
ing  the  ICP  magneto-hydrodynamic  equations  in  multi¬ 
dimensional  configurations.  Examples  are  given  in  the 
papers  by  Boulos,17  Mostaghimi,  Proulx  Bouos  and  co¬ 
workers, 18-23  Kolesnikov  and  co-workers, 24-30  Chen  and 
Pfender,31  van  den  Abeele  and  Degrez32,  and  more  re¬ 
cently  Panesi  et.  al.33  and  Morsli  and  Proulx.34  In  these 
works,  numerical  solutions  are  obtained  using  an  explicit 
coupling  approach,  by  solving  independently  the  flow  and 
the  electric  field  governing  equations,  and  updating  cou¬ 
pling  terms  after  each  iteration.  As  observed  by  van  den 
Abeele  and  Degrez,32  the  explicit  coupling  approach  pre¬ 
vents  the  use  of  Newton’s  method  during  the  first  it¬ 
erations  and  requires  to  resort  to  more  conservative  it¬ 
erative  techniques  ( e.g  Picard’s  method32)  at  the  begin¬ 
ning  of  the  simulation.  Convergence  issues  with  Newton’s 
method  may  also  arise  when  solving  the  flow  governing 
equations  in  time-dependent  form.  This  is  especially  true 
when  the  current  intensity  in  the  inductor  is  updated  to 
match  the  imposed  value  of  the  power  dissipated  in  the 
plasma.  The  cause  of  the  instability  is  most  probably 
due  to  the  lagged  update  of  the  Joule  heating  term  in  the 
energy  equation,  which  is  quadratic  in  the  electric  field 
amplitude. 17,31,32 

Most  of  the  simulations  available  in  the  litera¬ 
ture  assume  that  the  plasma  in  the  torch  is  in  LTE 
conditions.17-19,31-33,35  This  assumption  is  often  justified 
by  saying  that,  for  the  pressure  values  at  which  ICP  fa¬ 
cilities  are  operated  (e.g.  ss  104  Pa  and  above),  the  colli- 
sional  rate  among  the  gas  particles  are  sufficiently  large  to 
maintain  local  equilibrium.  A  second,  and  more  practical 
reason,  is  the  significant  stiffness  and  CPU  time  reduc¬ 
tion  compared  to  Non-LTE  (NLTE)  situations.  More¬ 
over,  in  LTE  conditions,  the  gas  thermodynamic  and 
transport  properties  are  only  function  of  two  independent 
state  variables36  (e.g.  pressure  and  temperature) .  Hence, 
they  can  be  easily  tabulated  to  further  reduce  the  com¬ 
putational  time.  Simulations  performed  by  Mostaghimi 
et  al.20,22  and  by  Zhang  et  al.37  (in  Argon  and  air  plas¬ 
mas,  respectively)  have  shown,  however,  that  the  use  of 
the  LTE  assumption  may  not  always  hold. 

An  accurate  modeling  of  NLTE  effects  in  ICP  plas¬ 
mas  can  be  achieved  by  means  of  State-to-State  (StS) 
models.38-50  These  treat  each  internal  energy  state  as  a 
separate  pseudo  species ,  thus  allowing  for  non-Boltzmann 
distributions.  Rate  coefficients  are  usually  obtained 
through  quantum  chemistry  calculations51-56  or  through 
phenomenological  models  providing  a  simplified  descrip¬ 
tion  of  the  kinetic  process  under  investigation.57,58  State- 
to-State  models  provide  a  superior  description  compared 
to  conventional  multi-temperature  models,  which  are 


based  on  Maxwell-Boltzmann  distributions.59-62  How¬ 
ever,  due  to  the  large  number  of  governing  equations  to 
be  solved,  their  application  to  multi-dimensional  prob¬ 
lems  can  become  computationally  demanding.63-67 

The  purpose  of  the  present  paper  is  development  of  a 
tightly  coupled  non-equilibrium  model  for  ICP  RF  plas¬ 
mas.  To  alleviate  the  possible  occurrence  of  numerical 
instabilities,  typical  of  an  explicit  coupling  approach,  the 
following  steps  are  taken:  (i)  the  flow  and  the  induced 
electric  field  governing  equations  are  solved  in  a  fully  cou¬ 
pled  fashion,  and  (ii)  steady-state  solutions  are  obtained 
by  means  of  a  time-marching  approach  (as  often  done  in 
CFD  applications68).  The  governing  equations  are  dis¬ 
cretized  in  space  by  using  the  Finite  Volume  method. 
Time-integration  is  then  performed  by  means  of  a  fully 
implicit  method.  As  it  is  shown  in  the  paper,  the  time- 
dependent  formulation  introduces  a  local  relaxation  in 
the  set  of  space-discretized  equations,  which  enhances 
convergence  significantly.  The  computational  ICP  frame¬ 
work  developed  in  this  work  allows  for  the  use  of  both 
LTE  and  NLTE  physico-chemical  models.  Non-LTE  ef¬ 
fects  are  described  based  on  a  hybrid  StS  model.  A  multi¬ 
temperature  (MT)  formulation  is  used  to  account  for 
thermal  non-equilibrium  between  translation  of  heavy- 
particles  and  vibration  of  molecules.  Excited  electronic 
states  of  atoms  are  instead  treated  as  separate  pseudo¬ 
species.  Free- electrons  are  assumed  Maxwellian  at  their 
own  temperature. 

The  paper  is  organized  as  follows.  Section  II  describes 
the  physical  model.  The  numerical  method  for  solving 
the  governing  equations  is  given  in  Sec.  III.  Computa¬ 
tional  results  are  presented  in  Sec.  IV.  Conclusions  are 
discussed  in  Sec.  V. 


II.  PHYSICAL  MODELING 

This  section  describes  the  physical  model  developed 
for  the  investigation  of  non-equilibrium  effects  in  ICP 
RF  plasmas.  The  non-equilibrium  model  for  the  ICP 
torch  is  built  based  on  the  torch  geometry  displayed  in 
Fig.  2.  To  make  the  problem  tractable,  the  following 
assumptions  are  introduced: 

(i)  Constant  pressure  and  no  macroscopic  streaming , 

(ii)  Charge  neutrality  and  no  displacement  current, 

(iii)  Steady-state  conditions  for  gas  quantities  (i.e. 

d()/dt  =  o), 

(iv)  No  gradients  along  the  axial  and  circumferential  di¬ 
rections  (i.e.  d()/dz  =  0,  d()/dc/)  =  0). 
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The  boundary  condition  at  the  torch  wall  is  obtained  as 
follows.  The  amplitudes  of  the  toroidal  electric  field  and 
the  axial  magnetic  field  are  linked  via:6 


1  drE 
r  dr 


—uuB, 


(9) 


FIG.  2.  Torch  geometry  and  adopted  reference  frame. 


A.  Electro-magnetic  field 


The  electromagnetic  field  inside  the  ICP  torch  is  de¬ 
scribed  by  the  Maxwell  equations: 


V  •  E  =  — ,  V  •  B  =  0,  (1) 

£o 

(9B  d~E 

V  x  E  =  — V  x  B  =  fj, qJ  +  (2) 


dt 


where  quantities  E  and  B  are  the  electric  and  magnetic 
fields,  respectively.  Quantity  pc  stands  for  the  charge 
density.  The  current  density  J  is  assumed  to  obey  Ohm’s 
law  J  =  creE,  with  ae  being  the  electrical  conductiv¬ 
ity.  Quantities  eo  and  po  are  the  vacuum  permittivity 
and  magnetic  permeability,  respectively.  The  applica¬ 
tion  of  the  simplifying  assumptions  just  introduced  to  the 
Maxwell  equations  (l)-(2)  leads  to  the  induction  equation 
for  the  induced  toroidal  electric  field: 


d  ( 1  drE 0  \  _  dE 0 

dr  \r  dr  )  ~  Mo<Te  dt  ' 


(3) 


Since  the  induced  eddy  currents  which  are  responsible 
for  the  heating  of  the  gas  are  induced  by  a  primary 
current  whose  intensity  varies  sinusoidally  in  time,  it 
seems  natural  to  seek  for  a  monochromatic  wave  solution, 
E^  =  Eexp(iwt),  where  w  =  2nf  (with  /  being  the  fre¬ 
quency  of  the  primary  current).  To  account  for  the  pos¬ 
sible  phase  difference  between  the  electric  and  magnetic 
fields,  the  amplitude  E  is  taken  complex,  E  =  Ele  +  iEim. 
The  substitution  of  E  exp(*wf)  in  Eq.  (3)  leads  to: 


0  x 


5rUem 

dt 


drFem 

dr 


—  ^Sem. 


(4) 


The  electromagnetic  (em)  conservative  variable,  flux  and 
source  term  vectors  are: 


Ucm  —  Eie  E\ ,  n  , 


F  = 

em 


Q 

uem 


dEre  dEh 


1  T 


dr  dr 
E  E- 

J-'re  7-,  Mm  ^ 

— 5-  +  WfloaeEim  —2 - Ldfl()(JeEr( 


1  T 


(5) 

(6) 

(7) 


Equation  (4)  must  be  supplemented  with  boundary  con¬ 
ditions  at  the  axis  (r  =  0)  and  at  the  torch  wall  (r  =  i?, 
with  R  being  the  torch  radius).  On  the  axis,  due  to  sym¬ 
metry,  both  components  of  the  electric  field  must  vanish: 


(8) 


where  the  amplitude  B  is  taken  complex.  Immediately 
outside  the  wall  the  magnetic  field  must  be  real  and,  since 
there  is  no  plasma  outside  the  tube,  its  value  can  only 
depend  on  the  ICP  operating  conditions  and  characteris¬ 
tics.  If  the  torch  is  long  enough,  the  magnetic  field  at  the 
torch  wall  can  be  approximated  with  the  expression  for 
an  infinite  solenoid,  B  =  po NIC,  where  quantities  N  and 
Jc  are  the  number  of  turns  per  unit-length  and  the  am¬ 
plitude  of  the  primary  current.  The  evaluation  of  Eq.  (9) 
at  the  torch  wall  and  the  use  of  the  relation  B  =  poNIc 
gives  the  wall  boundary  condition  for  the  induced  electric 
field: 


1  drEre  _  Q  1  drEim 
r  dr  1  r  dr 


ujpoNIc,  at  r  =  R.  (10) 


B.  Hydro-dynamics 

The  gas  contained  in  the  torch  is  made  of  electrons, 
atoms  and  molecules.  Charged  particles  comprise  elec¬ 
trons  and  positively  singly  charged  ions.  The  set  S 
stores  the  chemical  components,  and  the  heavy-particles 
are  stored  in  the  set  S h.  The  atomic  and  molecular 
components  are  stored  in  the  sets  <Sa  and  <Sm,  respec¬ 
tively.  The  previously  introduced  sets  satisfy  the  rela¬ 
tions  iSh  =  *Sa  U  <Sm  and  <S  =  {e-}  U  <Sh,  where  the 
symbol  e_  indicates  the  free-electrons.  The  electronic 
levels  of  the  heavy  components  are  stored  in  sets  if 
(with  s  £  «Sh)  and  are  treated  as  separate  pseudo  species 
based  on  a  StS  approach.69  The  notation  s*  is  used  to 
denote  the  i-th  electronic  level  of  the  heavy  component 
s  £  <5>h,  with  the  related  degeneracy  and  energy  being 
gf  and  Ef,  respectively.  A  multi-temperature  model  is 
instead  used  for  vibration  of  molecules  and  translation 
of  free-electrons  (with  the  related  temperatures  being  Tv 
and  Te,  respectively).60  Rotational  non-equilibrium  ef¬ 
fects  are  disregarded. 


Thermodynamic  properties 

The  gas  pressure  is  computed  as  p  =  pe+ph,  where  the 
symbol  ke  stands  for  Boltzmann’s  constant.  The  par¬ 
tial  pressures  of  free-electrons  and  heavy-particles  are, 
respectively,  pe  =  nekBTe  and  p^  =  nhkgT,  where  quan¬ 
tities  ne  and  Uh  denote,  respectively,  the  related  the  num¬ 
ber  densities  (with  n h  =  "}2sesh  ns  and  ns  =  nsJ- 

The  gas  total,  rotational,  vibrational  and  free-electron 
energy  densities  are: 

Pe=\p  +  per  +Pe''+E  E  ns '  +  A£s)  (n) 

seSh  iexf 


Eie  =  0,  E{m  =  0,  at  r  =  0. 
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per  =  ^  n^s(T)i  Pev  =  ns-^KTv),  pee  =  ^ Pe ■ 

s£<Sm  sG*Sm 

(12) 

Quantity  A.E1  stands  for  the  formation  energy  (per  par¬ 
ticle)  of  the  heavy  component  s  £  Sh-  The  average 
particle  rotational  and  vibrational  energies  (ETS  and  E J, 
s  £  Sm  Respectively)  are  computed,  respectively,  accord¬ 
ing  to  the  rigid-rotor  and  harmonic-oscillator  models70. 
Thermodynamic  data  used  in  this  work  are  taken  from 
Gurvich  tables71  (with  the  exception  of  the  spectroscopic 
data  for  the  electronic  levels  taken  from  Ref.  43). 


Chemical-kinetics 

The  NLTE  kinetic  mechanism  for  ICP  RF  plasmas  de¬ 
veloped  in  this  work  accounts  for: 

(i)  Excitation  by  electron  impact, 

(ii)  Ionization  by  electron  impact, 

(iii)  Dissociation  by  electron  impact, 

(iv)  Dissociation  by  heavy-particle  impact, 

(v)  Associative  ionization. 

The  endothermic  rate  coefficients  for  electron  induced 
processes  and  associative  ionization  reactions  are  taken 
from  the  ABBA  StS  model. 43-46  Those  for  dissociation  by 
heavy-particle  impact  are  taken  from  the  work  of  Park.60 
Reverse  rate  coefficients  are  obtained  based  on  micro¬ 
reversibility.  72,73 

The  mass  production  terms  for  free-electrons  and 
heavy-particles  are  computed  based  on  the  zeroth-order 
reaction  rate  theory.72,73  In  what  follows,  the  latter  quan¬ 
tities  are  indicated  with  the  notation  uje  and  los.  ,  respec¬ 
tively 

The  energy  transfer  terms  for  the  gas  vibrational  en¬ 
ergy  account  for  (i)  vibrational-translational  (vt)  en¬ 
ergy  exchange  in  molecule  heavy-particle  collisions,  (ii) 
vibrational-electron  (ve)  energy  exchange  in  molecule 
electron  collisions,  and  (iii)  the  creation/destruction  of 
vibrational  energy  in  chemical  reactions  (cv).  The  first 
two  energy  transfer  terms  (indicated  in  what  follows 
with  flvt  and  f lve,  respectively)  are  evaluated  based  on 
a  Landau- Teller  model,74  while  the  chemistry-vibration 
coupling  term  (flcv)  is  computed  by  using  the  non- 
preferential  dissociation  model  of  Candler.75  The  re¬ 
laxation  times  for  vt  energy  transfer  are  computed  by 
means  of  the  modified  formula  of  Millikan  and  White 
proposed  by  Park.60  The  energy  transfer  in  molecule- 
electron  inelastic  collisions  is  considered  only  for  N2.  The 
corresponding  relaxation  time  is  taken  from  the  work 
of  Bourdon.76  The  energy  transfer  terms  for  the  free- 
electron  gas  account  for  energy  exchange  undergone  by 
free-electrons  in  (i)  elastic  collisions  with  heavy-particles 
(f2ei),  (ii)  inelastic  electron  induced  excitation,  ioniza¬ 
tion  and  dissociation  processes  (flin)  and  (iii)  Joule  heat¬ 
ing  (Dj).  The  expressions  for  the  first  two  can  be 


found  in  Refs.  44-46.  The  (time-averaged)  Joule  heat¬ 
ing  source  term  is  obtained  by  averaging  over  a  pe¬ 
riod  the  instantaneous  Joule  heating  power  and  reads 
Dj=CTe(JE;r2e  +  E2n)/2.31-32 


Transport  properties  and  fluxes 

Transport  phenomena  are  treated  based  on  the  re¬ 
sults  of  the  Chapman-Enskog  method  for  the  Boltzmann 
equation77  under  the  assumption  that:  (i)  inelastic  and 
reactive  collisions  have  a  no  effect  on  the  transport  prop¬ 
erties  and  fluxes  and  (ii)  the  collision  cross-sections  for 
elastic  scattering  do  not  depend  on  the  internal  quantum 
states. 

The  translational  component  of  thermal  conductivity 
is  At  =  ]Cse.Sh  asX>"  where  the  mole  fractions  of  the 
heavy  components  are  Xs  =  risk^T/p  ( s  £  S h).  The  co¬ 
efficients  CTg  are  solution  of  the  linear  (symmetric)  trans¬ 
port  system  for  the  translational  thermal  conductivity 
(see,  for  instance,  Ref.  72  for  the  details).  The  con¬ 
tributions  of  the  gas  rotational  and  vibrational  degrees 
of  freedom  to  the  thermal  conductivity  (Ar  and  Av,  re¬ 
spectively)  are  taken  into  account  by  means  of  the  gener¬ 
alized  Eucken’s  correction.72  The  thermal  and  electrical 
conductivity  of  the  electron  gas  are:78,79 

75  /27rkBTe  XeA 22 

e  8  BV  me  A1e1Ae22-(A1e2)2’ 

_  3  e2  [2nk^  XeA 11 

CTe“  2kBVmeTeA00A1e1-(A10)9’ 

where  the  mole  fraction  of  free-electrons  is  Xe  = 
nckATe/p  and  e  =  1.602  x  10~19  C  is  the  electron  charge. 
Quantities  A lJe  denote  the  Devoto  collision  integrals.78 

The  mass  diffusion  fluxes  are  found  by  solving  the 
Stefan-Maxwell  equations  under  the  constraints  of  global 
mass  conservation  and  ambipolar  diffusion.'8-81  The  dif¬ 
fusion  driving  forces  include  only  mole  fraction  gradi¬ 
ents.  In  view  of  the  assumed  independence  of  the  elas¬ 
tic  collision  cross-section  on  the  internal  quantum  states, 
the  Stefan-Maxwell  equations  are  solved  for  the  diffusion 
fluxes  of  chemical  components  (Je  and  Js,  s  £  <Sh,  respec¬ 
tively).  The  mass  diffusion  fluxes  for  the  internal  levels 
( JSi ),  are  then  found  as  shown  in  Ref.  6.  The  gas  total, 
rotational,  vibrational  and  free-electron  heat  flux  are: 


(13) 

(14) 
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Governing  equations 


The  right-hand-side  residual  reads: 


The  governing  equations  for  the  gas  chemical  compo¬ 
sition  and  temperature  distribution  in  the  ICP  torch  are: 


drUg 

dt 


dr  Fg 

dr 


rSg. 


(19) 


The  gas  (g)  conservative  variable,  flux  and  source  term 
vectors  are: 


Ug  -  [  Pe  Psi  P ^  P^v  P&e  ]  5 

(20) 

Fg  =  [  Je  JSi  Q  Qv  Qe  ]  > 

(21) 

Sg  —  ^  hJq  J  5 

(22) 

Z  £  Xg  7  S  £  ,  AVltjll  —  ^vt  — I-  ^ve  — I-  ^cv 

and  Cle  = 

Hel  T  Hjn  +  • 

The  boundary  conditions  used  for  solving  Eq.  (19) 
are  a  symmetry  boundary  condition  at  the  axis  and  an 
isothermal  non-catalytic  boundary  condition  at  the  torch 
wall.  Following  the  work  of  Mostaghimi  et  al., 20,22  an  adi¬ 
abatic  wall  boundary  condition  is  used  for  the  vibrational 
and  free-electron  temperatures. 


Res j  =  rj+  iFj+i  -  rj_iFj_i  -  S j  ry  Ary,  (26) 

where  the  cell  volume  (length)  and  its  centroid  lo¬ 
cation  are  computed  as  Ary  =  ry+1/ 2  —  rj-i/2  and 
ry  =  (rj+ 1/2  +ri-i/2)/2,  respectively.  The  evaluation  of 
the  diffusive  flux  Fj+1/2  is  performed  by  approximating 
the  values  and  the  gradients  of  a  given  quantity  p  ( e.g . 
temperatures  and  electric  field  components)  by  means 
of  an  arithmetic  average  and  a  second  order  central  fi¬ 
nite  difference,  respectively.  To  facilitate  the  implemen¬ 
tation  of  the  constant  pressure  constraint,  the  solution 
update  is  performed  on  primitive  variables  (P)  consisting 
of  mass  fractions,  temperatures  and  electric  field  compo¬ 
nents,  P  —  (ye:ySi1  T ,  Tv,  -^e,  F/re,  Fim). 

dP 

FTj  ^  7'j Ary  =  —Res.,.  (27) 

The  transformation  matrix  T  can  be  obtained  from  the 
time- derivative  of  the  conservative  variables  (dXJ/dt)  by 
exploiting  the  global  continuity  equation,  dp/dt  =  0. 


III.  NUMERICAL  METHOD 


The  governing  equations  for  the  gas  and  the  electro¬ 
magnetic  fields  are  strongly  coupled  due  to  the  presence 
of  the  Joule  heating  term  in  Eq.  (22)  and  the  electrical 
conductivity  in  Eq.  (7).  This  suggests  to  adopt  a  fully 
coupled  approach  by  casting  Eqs.  (4)  and  (19): 


<9rrU 

dt 


dr  F 
dr 


=  rS, 


(23) 


where  the  conservative  variable,  flux  and  source  term 
vectors  are  now  U  =  (Ug,  Uem),  F  =  (Fg,  Fem)  and 
S  =  (Sg,  Sem),  respectively.  The  matrix  T  in  Eq.  (23) 
reads: 


A.  Temporal  discretization 

Equation  (27)  is  integrated  in  time  by  means  of  the 
backward  Euler  method:68 

<5P" 

rT^ryAry  =  -Res”+1,  (28) 

where  JP"  =  P”+1  —  P” .  The  local  time-step  At  is 
computed  based  on  the  von  Neumann  number  (£)  as 
At  =  £/(2 /?d),  where  quantity  pd  stands  for  the  spectral 
radius  of  the  diffusive  flux  Jacobian  matrix  <9F/9U.68  To 
advance  the  solution  from  the  time-level  n  to  the  time- 
level  n+1,  Eq.  (28)  is  linearized  around  the  time-level 
n.  The  outcome  of  the  linearization  is  a  block-tridiagonal 
algebraic  system  to  be  solved  at  each  time-step:68,82 


_  (  ^(ns+nt)  X  (ns+nt)  ^(ns-j-nt)x2 

\  02x(ns+nt)  02X2 

where  quantities  I  and  0  are,  respectively,  the  identity 
and  null  matrices.  Their  number  of  rows  and  columns 
are  indicated  by  the  first  and  second  lower-scripts,  re¬ 
spectively.  The  symbols  ns  and  nt  denote,  respectively, 
the  number  of  species  and  temperatures. 


Spatial  discretization 


L”  SP]^  +  C]  6P]  +  R”  <5P?+1  =  —Res”,  (29) 
where  the  left,  right  and  central  block  matrices  are: 

-A,-_i ,  (30) 

(31) 


L:j  = 

0  A  r 


2rj- 1/2 


j  T  Ary_i 
2  rj+ 1/2 


0-0 


T?  = _ _ A  , 

*0  A„  i  A  v 


A  ry  +  A?y+1 " 


Ci  = 


ft, 

/5S\  1 

<1 

_ 1 

;Ary  (Lj+Rj),  (32) 


The  application  of  the  Finite  Volume  method  to  Eq. 
(23)  leads  to  the  following  ODE  governing  the  time- 
evolution  of  the  conservative  variables  of  cell  j:68 

d\j  ■ 

F  q_i_  rj  Ary-  =  — Resj.  (25) 


where  quantity  A  stands  for  the  (primitive)  diffusive  flux 
Jacobian  matrix,  A  =  dF/dP.  The  block-tridiagonal 
system  (29)  is  solved  by  means  of  Thomas’  algorithm68 
and  the  solution  updated  at  the  next  time-level,  P"+1  = 
P”  +  <5P” .  This  process  is  continued  until  steady-state 
is  reached. 
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Notice  that  the  discretized  time  derivative  in  Eq.  (32) 
plays  the  role  of  a  relaxation  term.  Setting  this  term 
to  zero  ( i.e .  infinite  time-step)  is  equivalent  to  solve 
the  steady-state  form  of  the  governing  equations  (23) 
by  means  of  Newton’s  method.  Preliminary  calculations 
performed  in  LTE  conditions  indicated  that  this  strategy 
can  easily  lead  to  numerical  instability  problems  in  the 
initial  transient  of  the  simulation. 


IV.  APPLICATIONS 

The  gas  contained  in  torch  consists  of  molecular  nitro¬ 
gen  and  the  related  dissociation  and  ionization  products, 
S  =  {e_,  N,  N2,  Nj",  N+}.  Simulations  have  been  per¬ 
formed  by  means  of  the  ABBA  StS  model43-46  and  the 
MT  model  developed  by  Park.60 

The  torch  radius,  the  number  of  coils  per  unit-length 
and  the  wall  temperature  are  set  to  0.08  m,  50,  and  350  K, 
respectively.  The  current  intensity  of  the  primary  circuit 
is  found  from  the  solution  by  imposing  that  the  dissipated 
power  (per  unit-length)  in  the  plasma: 

P  =  2n  j  f2j  r  dr,  (33) 

J  o 

is  equal  to  a  fixed  value  Pq.  To  match  the  condition 
P  =  Pq  at  steady-state,  the  current  intensity  is  multi¬ 
plied  by  the  scaling  factor  7  =  y/ Pq/P  after  updating 
the  solution  at  the  new  time-step.  This  approach  was 
originally  introduced  by  Boulos17  in  the  1970’s,  and  it 
has  been  used  since  then  by  other  investigators. 18,31-33,37 
In  the  present  work,  no  scaling  is  applied  to  the  electric 
field  (as  opposed  to  explicit  coupling  methods32),  due  to 
the  use  of  a  fully  coupled  formulation. 

In  order  to  assess  the  influence  of  the  ICP  operating 
conditions  on  non-equilibrium  effects,  different  values  of 
pressure,  frequency  of  the  primary  circuit  and  dissipated 
power  have  been  adopted  (see  Table  I). 

TABLE  I.  Adopted  values  for  the  pressure,  frequency  of  the 
primary  circuit  and  dissipated  power  (per  unit-length)  in  the 
plasma. 


Quantity 

Units 

Values 

Pressure  (p) 

Pa 

3000,  5000  and  10  000 

Frequency  (/) 

MHz 

0.1,  0.5,  1  and  2.5 

Dissipated  power  (Pq) 

MW/m 

0.2,  0.3,  0.35  and  0.4 

A.  Performance  of  the  coupled  numerical  formulation 

The  coupled  numerical  method  developed  in  Sect.  Ill 
has  been  applied  over  a  wide  spectrum  of  operating  con¬ 
ditions  (given  in  Table  I)  in  both  LTE  and  NLTE  condi¬ 
tions.  In  all  cases  treated  in  this  work,  no  need  arose  for 
the  use  of  techniques  to  cope  with  numerical  instabilities 


such  as  under-relaxation  factors,  conservative  iterative 
strategies  ( e.g .  Picard’s  method)  or  choice  of  a  smart 
initial  guess  for  the  solution. 

To  demonstrate  in  practice  the  robustness  and  the 
effectiveness  of  the  proposed  computational  method,  a 
MT  NLTE  simulation  has  been  chosen  as  working  ex¬ 
ample.  The  operating  conditions  are:  p  =  10  000  Pa, 
/  =  0.5  MHz  and  Po  =  0.35  MW/m.  For  the  sake  of  sim¬ 
plicity,  an  isothermal  wall  boundary  conditions  has  been 
used  for  all  temperatures.  The  solution  has  been  initial¬ 
ized  with  a  uniform  equilibrium  distribution  at  7500  K, 
with  both  the  real  and  imaginary  electric  field  compo¬ 
nents  set  to  0.1  V/m.  The  initial  value  of  the  current 
intensity  was  250  A. 


FIG.  3.  Time-history  of  (a)  current  intensity  and  (b)  nor¬ 
malized  1/2  norm  of  the  relative  error  on  the  translational- 
rotational  temperature  for  the  MT  NLTE  model  ( p  = 
10  000  Pa,  /  =  0.5  MHz,  Po  =  0.35  MW/m;  isothermal  bound¬ 
ary  condition  for  all  temperatures). 

In  the  early  stages  of  the  numerical  simulation,  strong 
gradients  in  temperature  and  chemical  composition  form 
in  correspondence  of  the  wall,  due  to  the  isothermal 
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boundary  condition  imposed.  Despite  the  challenges  im¬ 
posed  by  the  problem,  the  initial  value  of  the  von  Neuman 
number  was  set  to  to  1  x  105,  and  increased  interactively 
(every  25  iterations)  up  to  5  x  108.  Such  large  values 
were  possible  due  to  the  adoption  of  a  fully  implicit  time- 
integration  method.  Figure  3  shows  the  time-history  of 
the  current  intensity  and  the  L2  norm  of  the  relative  error 
on  the  translational-rotational  temperature: 


errj,(r)  = 


\ 


3= 1 


(34) 


where  quantity  nc  denotes  the  number  of  cells,  and 
ST™  =  T"+1  —  T".  A  converged  solution  is  achieved 
in  less  than  100  iterations  (see  Fig.  3(b)),  with  a  reduc¬ 
tion  of  more  than  nine  orders  of  magnitude  in  the  relative 
error  for  the  temperature.  In  practice,  the  solution  is  al¬ 
ready  converged  after  80  iterations,  as  can  be  observed 
from  the  current  intensity  time-history  (see  Fig.  3(a)). 

Figure  4  shows  the  temperature  distribution  along 
the  torch  radius.  Due  to  the  efficient  energy  exchange 
in  N2-e“  interactions,60,76,83  the  vibrational  and  free- 
electron  temperature  profiles  are  essentially  indistin¬ 
guishable.  This  feature  has  been  observed  in  all  NLTE 
simulations  performed.  This  is  the  reason  for  the  use  of 
the  notation  Tve  in  Fig.  4  (and  also  in  what  follows). 


FIG.  4.  MT  NLTE  temperature  distribution  ( p  =  10  000  Pa, 
/  =  0.5  MHz,  Po  =  0.35  MW/m;  unbroken  line  T,  dashed-line 
Tve;  isothermal  boundary  condition  for  all  temperatures). 


B.  Assessment  of  non-equilibrium  effects 

Before  discussing  in  detail  NLTE  effects,  LTE  simula¬ 
tions  have  been  performed  and  compared  with  the  MT 
NLTE  results  to  assess  the  extent  of  the  departure  from 
equilibrium.  The  pressure  is  set  to  lOOOOPa  as,  for  this 
relatively  high  value,  LTE  conditions  are  often  assumed.6 
The  frequency  and  the  dissipated  power  are  0.5  MHz  and 
0.35  MW/m,  respectively. 


r(m) 

(b) 


FIG.  5.  Comparison  between  the  LTE  and  the  MT  NLTE 
temperature  (a)  and  Joule  heating  (b)  distributions  ( p  = 
lOOOOPa,  /  =  0.5 MHz,  Po  =  0.35 MW/m;  unbroken  line 
LTE,  dashed  line  NLTE). 


Figure  5  compares  the  LTE  and  NLTE  translational- 
rotational  temperature  and  Joule  heating  distributions. 
Close  to  the  wall,  the  curvature  of  the  LTE  tempera¬ 
ture  distribution  changes  sign.  This  is  a  consequence 
of  the  non-monotone  behavior  of  the  equilibrium  total 
thermal  conductivity  of  the  working  gas  (nitrogen).  The 
same  trend  is  also  found  for  the  MT  solution,  though  less 
pronounced.  The  MT  model  predicts  that  the  gas  is  in 
thermal  equilibrium  close  to  the  axis  (see  Fig.  4).  In 
both  the  LTE  and  NLTE  simulations,  the  temperature 
is  maximum  on  the  axis,  due  to  the  absence  of  radiative 
losses  in  the  plasma.9,10  Overall,  the  LTE  solution  pre¬ 
dicts  higher  temperature  values,  with  the  difference  being 
maximum  on  the  axis.  This  is  a  general  trend  observed  in 
all  the  simulations  performed  in  this  work  (and  also  in  the 
multi-dimensional  results  obtained  by  other  investigators 
in  Refs.  6  and  37).  In  NLTE  conditions,  the  temperature 
is  lower  because  the  plasma  is  heated  over  a  wider  region 
compared  to  LTE  conditions.  This  is  confirmed  by  the 


Joule  heating  distribution  shown  in  Fig.  5(c).  Analogous 
conclusions  can  be  drawn  when  comparing  with  the  StS 
NLTE  model  adopted  in  this  work. 

The  results  in  Fig.  5  demonstrate  that,  even  at  rela¬ 
tively  high  pressures,  the  LTE  assumption  can  lead  to  a 
severe  overestimation  of  the  gas  temperature  and,  in  gen¬ 
eral,  to  an  inaccurate  prediction  of  the  thermo-chemical 
state  of  the  gas.  This  is  further  confirmed  by  the  com¬ 
parison  for  the  mole  fractions  given  in  Fig.  6.  It  is 
worth  recalling  that  in  the  present  work,  the  effects  of 
macroscopic  flows  (which  enhance  non-equilibrium)  are 
neglected. 


FIG.  6.  Comparison  between  the  LTE  and  the  MT  NLTE 
mole  fraction  distributions  (p  =  10  000  Pa,  /  =  0.5  MHz,  Po  = 
0.35  MW/m;  unbroken  line  LTE,  dashed  line  NLTE).  Due  to 
the  low  concentration  of  N^  (not  shown),  the  mole  fractions 
of  e“  and  N+  are  essentially  indistinguishable  in  both  LTE 
and  NLTE  conditions. 


C.  Influence  of  frequency  and  power  on  non-equilibrium 

In  order  to  assess  the  influence  of  operating  conditions 
on  non-equilibrium  phenomena,  it  was  decided  to  per¬ 
form  a  parametric  study  on  the  frequency  of  the  primary 
circuit  and  the  power  dissipated  in  the  plasma.  Figures  7- 
8  show  the  results  of  this  investigation  for  the  MT  NLTE 
model  at  p  =  5000  Pa.  Increasing  the  frequency  (Fig.  7) 
has  the  effect  of  narrowing  the  extent  of  the  skin-depth , 
which  is  the  zone  over  which  most  of  the  power  is  dissi¬ 
pated  by  Ohmic  heating.  This  can  be  seen  from  the  re¬ 
lation  for  the  skin-depth14,15 ,  <5  =  (cre7'‘Mo/)_1,/2-  The  re¬ 
duction  of  the  skin-depth  induces  sharper  gradients  close 
to  the  wall,  thereby  enhancing  non-equilibrium.  These 
findings  are  in  accordance  with  the  observations  reported 
by  Mostaghimi  et  al.22  for  NLTE  Argon  plasmas.  The  in¬ 
crease  of  the  dissipated  power  (Fig.  8)  has  an  opposite 
effect.  As  can  be  observed  from  the  results,  higher  power 
levels  favor  the  establishment  of  thermal  equilibrium  con¬ 
ditions. 


FIG.  7.  MT  NLTE  temperature  distributions  for  different 
values  of  the  frequency  of  the  primary  circuit:  unbroken  line 
/  =  0.1  MHz,  dashed  line  /  =  0.5  MHz,  dotted-dashed  line 
/  =  1  MHz,  dotted  line  /  =  2.5  MHz  ( p  =  5000  Pa,  Po  = 
0.35  MW/m;  unbroken  line  T,  line  with  symbols  Tve). 


r  (m) 

FIG.  8.  MT  NLTE  temperature  distributions  for  differ¬ 
ent  values  of  the  dissipated  power  (per  unit-length)  in  the 
plasma:  unbroken  line  Po  =  0.2  MW/m,  dashed  line  Po  = 
0.3MW/m,  dotted-dashed  line  Po  =  0.35 MW/m,  dotted  line 
Po  =  0.4  MW/m  ( p  =  5000  Pa,  /  =  0.5  MHz;  unbroken  line 
T,  line  with  symbols  Tve). 


D.  Comparison  between  the  StS  and  the  MT  predictions 

This  section  compares  the  predictions  obtained  by  the 
StS  and  MT  NLTE  models.  Figure  9  shows  the  compar¬ 
ison  in  terms  of  temperatures  and  N  mole  fraction.  For 
both  the  StS  and  MT  solutions,  decreasing  the  pressure 
has  the  effect  of  enhancing  thermal  non-equilibrium.  In 
the  case  of  the  MT  model,  the  translational-rotational 
temperature  is  maximum  on  the  axis  and  decreases 
monotonically  when  approaching  the  wall. 

On  the  other  hand  the  free-electron  temperature 
increases,  reaches  a  maximum  and  then  decreases 
until  it  reaches  the  value  determined  from  the  adiabatic 
bound- 
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FIG.  9.  Comparison  between  the  translational-rotational  temperature  (left)  and  the  N  mole  fraction  (right)  distributions 
predicted  by  the  StS  and  the  MT  NLTE  models  at  different  pressures:  (a)-(b)  p  =  10  000  Pa,  (c)-(d)  p  =  5000  Pa,  (e)-(f) 
p  =  3000 Pa  (/  =  0.5 MHz,  Po  =  0.35 MW/m;  in  (a),  (c)  and  (e)  unbroken  line  T  StS,  dashed  line  Tve  StS,  dotted-dashed  line 
T  MT,  dotted  line  Tve  MT;  in  (b),  (d)  and  (f)  unbroken  line  StS,  dashed  line  MT). 


ary  condition.  This  behavior  is  due  to  the  balance  be¬ 
tween  the  Joule  heating  (which  heats  up  the  electron  gas) 
and  the  energy  loss  in  elastic  and  inelastic  collisions,  and 


chemical  reactions.  The  peak  location  of  the  free-electron 
temperature  moves  towards  the  wall  when  decreasing  the 
pressure.  This  is  a  consequence  of  the  Joule  heating  dis- 
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tribution  (not  shown  in  Fig.  9)  which  becomes  sharper 
and  clustered  to  the  wall  at  lower  pressures.  It  is  worth 
noticing  that,  in  the  StS  simulation,  (i)  the  translational- 
rotational  temperature  no  longer  exhibits  a  monotone 
behavior  and  (ii)  the  axis  values  of  both  temperatures 
are  systematically  lower  than  those  predicted  by  the  MT 
model.  The  differences  observed  between  the  StS  and 
MT  temperature  distribution  have  an  effect,  as  it  should 
be  expected,  on  the  chemical  composition  (see  Figs.  9(b), 
9(d)  and  9(f)). 

Figure  10  shows  the  normalized  population  of  the  elec¬ 
tronic  levels  of  N  on  the  torch  axis  (circles) ,  in  the  mid¬ 
point  of  the  torch  (squares)  and  at  the  wall  (triangles) 
at  different  pressures.  The  population  exhibit  significant 
distortions  from  a  Boltzmann  shape  only  close  to  the  wall 
(where  recombination  occurs).  Deviations  from  a  Boltz¬ 
mann  distributions  are  more  significant  when  increasing 
the  pressure  due  to  higher  recombination. 

V.  CONCLUSIONS 

A  tightly  coupled  magneto-hydrodynamic  solver  for 
the  study  of  the  weakly  ionized  plasmas  found  in  RF 
discharges  has  been  developed.  A  hierarchy  of  thermo¬ 
physical  models  have  been  added  to  the  solver  to  model 
the  non-equilibrium  effects  in  atomic  and  molecular  plas¬ 
mas.  These  include  LTE,  multi-temperature  and  the 
more  sophisticated  State-to-State  models.  The  govern¬ 
ing  equations  for  the  flow  and  electromagnetic  fields 
have  been  written  as  a  system  of  coupled  time-dependent 
conservation-laws.  Steady-state  solutions  have  been  ob¬ 
tained  by  means  of  an  implicit  Finite  Volume  Method. 

Results  obtained  by  using  a  multi-temperature  and 
State-to-State  models  have  shown  that  the  LTE  assump¬ 
tion  does  not  hold  and  that  its  use  can  lead  to  a  wrong 
prediction  of  the  thermo-chemical  state  of  the  gas.  The 
analysis  of  the  temperature  distribution  in  the  torch  indi¬ 
cated  that  non-equilibrium  plays  an  important  role  close 
to  the  walls,  due  to  the  combined  effects  of  Ohmic  heat¬ 
ing,  and  chemical  composition  and  temperature  gradi¬ 
ents.  The  accurate  study  of  the  population  of  excited 
electronic  states  has  shown  that,  in  view  of  the  absence 
of  a  macroscopic  gas  flow,  non-Boltzmann  distributions 
are  limited  to  a  narrow  region  close  to  the  torch  wall. 
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